MRG/ENH: apply_inverse_cov #6262

dengemann · 2019-05-02T18:50:44Z

Here we go: a candidate for a future apply_inverse_cov function that gets us source power by projecting the data covariance.

Feel free to edit!

cc @larsoner @britta-wstnr @wmvanvliet

lgtm-com · 2019-05-02T19:34:53Z

This pull request introduces 2 alerts when merging 9eced1e into 3cc76b9 - view on LGTM.com

new alerts:

2 for Unused import

Comment posted by LGTM.com

codecov · 2019-06-13T09:44:01Z

Codecov Report

Merging #6262 into master will decrease coverage by 4.42%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master    #6262      +/-   ##
==========================================
- Coverage   89.99%   85.57%   -4.43%     
==========================================
  Files         453      453              
  Lines       82020    82100      +80     
  Branches    13000    13002       +2     
==========================================
- Hits        73812    70254    -3558     
- Misses       5381     9177    +3796     
+ Partials     2827     2669     -158

wmvanvliet · 2019-06-13T09:47:16Z

Cool stuff @dengemann!

I've looked into how to combine the XYZ and your original was correct, except that when doing K @ cov @ K.T, the values are already squared, hence the combine_xyz should not square them again. I ended replacing the call to combine_xyz with just sol = sol[0::3] + sol[1::3] + sol[2::3].

Another change I made to the example is to also compute the power for the baseline period. Then, we can plot the power relative to baseline, which looks much cleaner.

bloyl · 2019-11-29T16:42:47Z

@larsoner @britta-wstnr @wmvanvliet @dengemann
Apart from restructuring and testing what needs to be done here?

dengemann · 2019-11-29T16:44:19Z

just that

…

On 29 Nov 2019, at 17:42, Luke Bloy ***@***.***> wrote: @larsoner @britta-wstnr @wmvanvliet @dengemann Apart from restructuring and testing what needs to be done here? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

bloyl · 2019-12-30T15:14:08Z

@larsoner @britta-wstnr @wmvanvliet @dengemann
Any opinions on testing? I could add a simple regression test or is there a 3rd party method/implementation we'd like to test against.

dengemann · 2019-12-30T15:23:01Z

A good one is to source localize with MNE and compute the covariance of the sources and retain the diagonal. The result is identical if all is done correctly. It may depend however on the solver settings. I can tell you though that it works with this code. You get the idea?

…

On 30 Dec 2019, at 16:14, Luke Bloy ***@***.***> wrote: @larsoner @britta-wstnr @wmvanvliet @dengemann Any opinions on testing? I could add a simple regression test or is there a 3rd party method/implementation we'd like to test against. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

bloyl · 2020-01-02T17:47:54Z

@larsoner @britta-wstnr @wmvanvliet @dengemann
This is ready for a code review. I can't change the title to [MRG] though.

A couple of potential issues remain and need to be discussed.

I'm only testing with pick_ori='normal' since that has a clear expected relationship to what is returned using apply_inverse_raw. I feel like that gets the coverage and the example provides shows good agreement between methods (LCMV and DICS) so I'm fine with an abbreviated test.
Using method='eLORETA' throws an exception in _assemble_kernel

mne-python/mne/minimum_norm/inverse.py

Lines 682 to 683 in 241e119

if method != "MNE":

noise_norm = noise_norm[src_sel]

which is strange since eLoreta seems to be tested in other parts of the code base. I don't know if this is a bug in _assemble_kernel or if eLoreta should not be an option for this method. Someone with knowledge of eLoreta should decide if it should be supported here.
The unit label for the colorbar in the topoplots of one of the example figures are overlaping with tick labels. https://17357-1301584-gh.circle-artifacts.com/0/dev/auto_examples/inverse/plot_mne_cov_power.html#sphx-glr-auto-examples-inverse-plot-mne-cov-power-py. Personally, I would remove the unit label but if there is another suggestion that's fine.
Someone with more familiarity with the MNE inverse filters should check that the docstring (mostly copied from apply_inverse) is correct for apply_inverse_cov.
Updating Whats_new is done. I'll put it in the next Push :)

mne/minimum_norm/tests/test_inverse.py

examples/inverse/plot_mne_cov_power.py

bloyl · 2020-02-16T15:35:18Z

This is ready for a final review.

agramfort

@bloyl if you could add more explanations to both examples eg by making sure you have titles on stc plot always and that you have a short explanation to explain what this figures show. Thx

doc/changes/latest.inc

examples/inverse/plot_evoked_ers_source_power.py

mne/minimum_norm/inverse.py

wmvanvliet · 2020-02-17T12:39:47Z

mne/minimum_norm/inverse.py

+    # always combine vector components
+    if is_free_ori:
+        logger.info('    Combining the current components...')
+        sol = sol[0::3] + sol[1::3] + sol[2::3]


So what happens when the user specifies pick_ori='vector'?

I guess I don't understand your question. If pick_ori='vector' then _is_free_ori is True and you combine the outputs from the inverse.

What happens with fixed orientation inverse operators is unfortunately not tested since you can only use pick_ori=None with those and I don't have a comparable call to apply_inverse_raw that will facilitate testing. Presumably, this would be the same as using pick_ori='normal' and a suitably configured inverse.

But pick_ori='vector' should produce a VectorSourceEstimate. We explicitly not support it in the apply_lcmv_cov and apply_dics_cov functions, as it should return vectors pointing in the direction of maximum power and we haven't gotten around to implementing that yet.

mne/minimum_norm/tests/test_inverse.py

bloyl · 2020-02-19T12:56:18Z

Im fine with disallowing it, my personal position is all the apply_inverse* methods should return vectors and the pick_ori functionality put in the source estimate objects. But that's another discussion. I will say that it would have been helpful to have this api conversation 3 months ago when I took over the pr and asked what restructuring needed to be done.

…

On Wed, Feb 19, 2020, 4:33 AM Marijn van Vliet ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In mne/minimum_norm/inverse.py <#6262 (comment)>: > + logger.info(' Picked %d channels from the data' % len(sel)) + logger.info(' Computing inverse...') + + K, noise_norm, vertno, source_nn = _assemble_kernel(inv, label, method, + pick_ori) + + # apply imaging kernel + sol = np.einsum('ij,ij->i', K, (cov.data[sel][:, sel] @ K.T).T)[:, None] + + is_free_ori = (inverse_operator['source_ori'] == + FIFF.FIFFV_MNE_FREE_ORI and pick_ori != 'normal') + + # always combine vector components + if is_free_ori: + logger.info(' Combining the current components...') + sol = sol[0::3] + sol[1::3] + sol[2::3] But pick_ori='vector' should produce a VectorSourceEstimate. We explicitly not support it in the apply_lcmv_cov and apply_dics_cov functions, as it should return vectors pointing in the direction of maximum power and we haven't gotten around to implementing that yet. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#6262?email_source=notifications&email_token=ABKTXHIA52HP5IWZTNGZZMDRDT4IPA5CNFSM4HKEEJT2YY3PNVWWK3TUL52HS4DFWFIHK3DMKJSXC5LFON2FEZLWNFSXPKTDN5WW2ZLOORPWSZGOCWB4XYA#discussion_r381172467>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABKTXHN77RQIGSLBL3MIKGTRDT4IPANCNFSM4HKEEJTQ> .

dengemann · 2020-02-19T14:37:54Z

Let me know once I should take a look / help.

larsoner · 2020-02-19T14:38:43Z

I will say that it would have been helpful to have this api conversation 3 months ago when I took over the pr and asked what restructuring needed to be done.

If you want me to look into the pick_ori='vector' changes I can give it a shot

bloyl · 2020-02-19T15:59:57Z

The main thing that needs to be decided is the api. @wmvanvliet has reasonably raised 2 questions:

Should we have a dB : bool transform data to decibels parameter. This was in the original PR and matches some(?) psd code. I prefer not to have data in decibels so I'm fine with either keeping the parameter or never returning decibels. I would say that whatever is decided it should be propagated to the other apply_*_cov methods.
What should we do with pick_ori=vector? I haven't thought too deeply about it but I feel like none of the apply_*_cov methods should return a VectorSourceEstimate since the whole k * cov * K.t formulation comes about from taking the inner product of source waveforms and thus collapsing orientation.
That being said, I believe that pick_ori=vector is functionally equivalent to pick_ori=None for free orientation inverses. So I have no issue with removing it as an allowed option.

Once those are decided they are minimal changes needed to the method, but the test will need to be reworked depending on the decisions.

dengemann · 2020-02-19T16:38:45Z

The main thing that needs to be decided is the api. @wmvanvliet has reasonably raised 2 questions:

I think 1 could be ok, although I'm not feeling very strongly about it. 2 is also a point that is somewhat confusing to me, I'd rather not support it at this point. Think that we have a power estimate here!

larsoner · 2020-02-19T16:49:23Z

I think we should get rid of the dB transformation and let users do it afterward if they want to.
The easiest path forward is probably to disallow pick_ori='vector' and we can add support for it in a subsequent PR.

I agree with @wmvanvliet that what you get out of pick_ori='vector' should be a VectorSourceEstimate, i.e. the three orientations should not be combined, whereas with pick_ori=None they should be combined (as they are everywhere else in MNE).

larsoner · 2020-03-02T19:52:55Z

Rebased, added support for pick_ori='vector', and simplified the code a bit by routing through apply_inverse itself (DRY).

The only remaining issue I see is the tolerance on the cov calc is not great, but that I think will be solved by #7369, let's try to get that in first, rebase this, and try again.

larsoner · 2020-03-03T12:54:56Z

@wmvanvliet feel free to take a look now, other than the rtol issue that should hopefully be fixed by #7369 we should be good to go.

bloyl · 2020-03-03T16:13:47Z

mne/minimum_norm/inverse.py

+    sol = cov.data[sel][:, sel] @ K.T
+    sol = np.sum(K * sol.T, axis=1, keepdims=True)


I feel like this is not equivalent to
np.diag(K @ cov.data[sel][:, sel] @ K.T)
Which i think is what we want. Doesn't the current code include off-diagonal elements for the src covariance matrix?

I feel like this is not equivalent to ...

It should be equivalent, though admittedly it is confusing because it looks like it might be computing something different. The diagonal entries of A @ B are given by (A * B.T).sum(axis=1), i.e., the sum of the first-row-of-A-times-first-col-of-B, sum of second-row-of-A-times-second-col-of-B, etc. Just to check it numerically, this code works fine locally if I add it:

sol = np.sum(K * sol.T, axis=1, keepdims=True) sol2 = np.diag(K @ cov.data[sel][:, sel] @ K.T)[:, np.newaxis] sol3 = np.sum((K @ cov.data[sel][:, sel]) * K, axis=1, keepdims=True) np.testing.assert_allclose(sol, sol2) np.testing.assert_allclose(sol, sol3)

Doesn't the current code include off-diagonal elements for the src covariance matrix?

Nope, just the diagonal. The end result is has shape (n_src, 1).

larsoner · 2020-03-04T00:32:00Z

... based on what @wmvanvliet said above:

But pick_ori='vector' should produce a VectorSourceEstimate. We explicitly not support it in the apply_lcmv_cov and apply_dics_cov functions, as it should return vectors pointing in the direction of maximum power and we haven't gotten around to implementing that yet.

I think it makes the most sense just not to allow pick_ori='vector', at least for now. We can always add support for it later.

Good to go for you in that case @bloyl ?

wmvanvliet · 2020-03-04T06:39:14Z

Going from the CSD to the max-power orientation is computing the largest eigenvector of the matrix. You can see how it is done in the beamformer code when pick_ori='max-power'. I'm fine with just disabling it for now. We can get to it later.

britta-wstnr

I am late to the party here, apologies. I added some small comments regarding the examples.

Regarding the almost philosophical discussion about vector solutions and power estimates: to me it makes sense to have one power estimate per source orientation in that case and I have used those before (mostly for checking data/forward models). That's my 2 cents.

examples/inverse/plot_evoked_ers_source_power.py

examples/inverse/plot_mne_cov_power.py

examples/inverse/plot_evoked_ers_source_power.py

examples/inverse/plot_mne_cov_power.py

mne/minimum_norm/inverse.py

mne/minimum_norm/tests/test_inverse.py

agramfort · 2020-03-04T12:10:43Z

then let's do: Here we will use dSPM as well as DICS and LCMV beamformers.

…

agramfort · 2020-03-04T14:35:04Z

ok !

…

larsoner · 2020-03-04T15:18:15Z

doc/glossary.rst

+    LCMV beamformer
+        Linearly constrained minimum variance beamformer, which attempts to
+        estimate activity for a given source while suppressing cross-talk from
+        other regions, see :func:`mne.beamformer.make_lcmv`.


@britta-wstnr @agramfort thanks for the comments, pushed a commit to address them. In the example I linked to :term:s instead of giving the definitions there, which caused me to add LCMV and DICS to our glossary.rst. @britta-wstnr can you just check to see if this is okay? If it's not, let me know how to modify it.

@britta-wstnr I'll merge this once CIs come back happy, but feel free to open a PR to correct the wording

@larsoner, thanks it is on my list to look at the glossary also with the beamformer tutorial I am writing!

FYI @britta-wstnr if you can get to the tutorial in the next couple of weeks it can make it into 0.20!

I will do all I can! It has been on my to do list for way too long already anyway 🙈 ... the workload has been high the last couple months.

* mne power example * Fix combining xyz, compute power relative to baseline * mv method to min_norm folder * Fix Doc strings - us dB instead of log * Add second comparitive example * Better description * Ref other example * WIP: failing test * Test passes of ori-normal fails for others * fix example cleanup test * Fix testing * rebase * disallow pick_ori=vector and rework test * restructure test * FIX: Marijn comments * ENH: Simplify by routing through Evoked * FIX: Fix examples * FIX: decim * DOC: Comment [ci skip] * FIX: Remove cruft * FIX: Fix after rebase * FIX: Remove vector * DOC: Better docs * DOC: Document another [ci skip] * FIX: Remove workaround Co-authored-by: Marijn van Vliet <[email protected]> Co-authored-by: Luke Bloy <[email protected]> Co-authored-by: Eric Larson <[email protected]>

dengemann mentioned this pull request May 2, 2019

MRG: Power estimate for LCMV: apply_lcmv_cov #6227

Merged

bloyl force-pushed the apply_inverse_cov branch from e403b04 to 78f624e Compare December 28, 2019 20:08

agramfort reviewed Jan 4, 2020

View reviewed changes

mne/minimum_norm/tests/test_inverse.py Outdated Show resolved Hide resolved

examples/inverse/plot_mne_cov_power.py Show resolved Hide resolved

bloyl force-pushed the apply_inverse_cov branch 2 times, most recently from 29418d0 to 330b8d3 Compare February 16, 2020 13:31

bloyl changed the title ~~WIP/ENH: apply_inverse_cov~~ MRG/ENH: apply_inverse_cov Feb 16, 2020

agramfort reviewed Feb 16, 2020

View reviewed changes

doc/changes/latest.inc Show resolved Hide resolved

wmvanvliet requested changes Feb 17, 2020

View reviewed changes

bloyl force-pushed the apply_inverse_cov branch from 330b8d3 to b51d7c3 Compare March 2, 2020 16:45

larsoner force-pushed the apply_inverse_cov branch from dcaf730 to a1ce007 Compare March 2, 2020 19:50

larsoner force-pushed the apply_inverse_cov branch from 21a4aba to 1654b78 Compare March 2, 2020 20:20

bloyl reviewed Mar 3, 2020

View reviewed changes

larsoner force-pushed the apply_inverse_cov branch from 04acb81 to a80c7c3 Compare March 3, 2020 16:35

bloyl and others added 10 commits March 3, 2020 18:10

rebase

d1f0e23

disallow pick_ori=vector and rework test

7fa8564

restructure test

e68cc43

FIX: Marijn comments

ccea55f

ENH: Simplify by routing through Evoked

8f57d99

FIX: Fix examples

bb7342b

FIX: decim

a099541

DOC: Comment [ci skip]

8239f18

FIX: Remove cruft

f23ce15

FIX: Fix after rebase

2ef5eaf

larsoner force-pushed the apply_inverse_cov branch from 00d2b84 to 2ef5eaf Compare March 3, 2020 23:14

FIX: Remove vector

ffee32a

wmvanvliet approved these changes Mar 4, 2020

View reviewed changes

britta-wstnr reviewed Mar 4, 2020

View reviewed changes

examples/inverse/plot_evoked_ers_source_power.py Outdated Show resolved Hide resolved

examples/inverse/plot_evoked_ers_source_power.py Outdated Show resolved Hide resolved

examples/inverse/plot_mne_cov_power.py Show resolved Hide resolved

agramfort reviewed Mar 4, 2020

View reviewed changes

larsoner mentioned this pull request Mar 4, 2020

ENH: Add cov.plot_topomap(info) #7384

Closed

larsoner added this to the 0.20 milestone Mar 4, 2020

DOC: Better docs

c15dfbe

larsoner reviewed Mar 4, 2020

View reviewed changes

larsoner mentioned this pull request Mar 4, 2020

WIP: Generate subject fiducials #6693

Closed

agramfort approved these changes Mar 4, 2020

View reviewed changes

larsoner added 2 commits March 4, 2020 10:46

DOC: Document another [ci skip]

3152e4e

FIX: Remove workaround

aca1135

larsoner merged commit c086e76 into mne-tools:master Mar 4, 2020

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

MRG/ENH: apply_inverse_cov #6262

MRG/ENH: apply_inverse_cov #6262

dengemann commented May 2, 2019 •

edited

Loading

lgtm-com bot commented May 2, 2019

codecov bot commented Jun 13, 2019 •

edited

Loading

wmvanvliet commented Jun 13, 2019 •

edited

Loading

bloyl commented Nov 29, 2019

dengemann commented Nov 29, 2019 via email

bloyl commented Dec 30, 2019

dengemann commented Dec 30, 2019 via email

bloyl commented Jan 2, 2020

bloyl commented Feb 16, 2020

agramfort left a comment

wmvanvliet Feb 17, 2020

bloyl Feb 18, 2020

wmvanvliet Feb 19, 2020

bloyl commented Feb 19, 2020 via email

dengemann commented Feb 19, 2020

larsoner commented Feb 19, 2020

bloyl commented Feb 19, 2020

dengemann commented Feb 19, 2020

larsoner commented Feb 19, 2020

larsoner commented Mar 2, 2020

larsoner commented Mar 3, 2020

bloyl Mar 3, 2020

larsoner Mar 3, 2020

larsoner commented Mar 4, 2020

wmvanvliet commented Mar 4, 2020

britta-wstnr left a comment

agramfort commented Mar 4, 2020 via email

agramfort commented Mar 4, 2020 via email

larsoner Mar 4, 2020

larsoner Mar 4, 2020

britta-wstnr Mar 5, 2020

larsoner Mar 5, 2020

britta-wstnr Mar 5, 2020

		sol = cov.data[sel][:, sel] @ K.T
		sol = np.sum(K * sol.T, axis=1, keepdims=True)

MRG/ENH: apply_inverse_cov #6262

MRG/ENH: apply_inverse_cov #6262

Conversation

dengemann commented May 2, 2019 • edited Loading

lgtm-com bot commented May 2, 2019

codecov bot commented Jun 13, 2019 • edited Loading

Codecov Report

wmvanvliet commented Jun 13, 2019 • edited Loading

bloyl commented Nov 29, 2019

dengemann commented Nov 29, 2019 via email

bloyl commented Dec 30, 2019

dengemann commented Dec 30, 2019 via email

bloyl commented Jan 2, 2020

bloyl commented Feb 16, 2020

agramfort left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bloyl commented Feb 19, 2020 via email

dengemann commented Feb 19, 2020

larsoner commented Feb 19, 2020

bloyl commented Feb 19, 2020

dengemann commented Feb 19, 2020

larsoner commented Feb 19, 2020

larsoner commented Mar 2, 2020

larsoner commented Mar 3, 2020

Choose a reason for hiding this comment

Choose a reason for hiding this comment

larsoner commented Mar 4, 2020

wmvanvliet commented Mar 4, 2020

britta-wstnr left a comment

Choose a reason for hiding this comment

agramfort commented Mar 4, 2020 via email

agramfort commented Mar 4, 2020 via email

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dengemann commented May 2, 2019 •

edited

Loading

codecov bot commented Jun 13, 2019 •

edited

Loading

wmvanvliet commented Jun 13, 2019 •

edited

Loading